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Abstract. We study soliton solutions to a nonlinear Schrodinger equation with a satu- 
rated nonlinearity. Such nonlinearities are known to possess minimal mass soliton solu- 
tions. We consider a small perturbation of a minimal mass soliton, and identify a system 
of ODEs similar to those from [16], which model the behavior of the perturbation for short 
times. We then provide numerical evidence that under this system of ODEs there are two 
possible dynamical outcomes, which is in accord with the conclusions of [31]. For initial 
QQ data which supports a soliton structure, a generic initial perturbation oscillates around 

I the stable family of solitons. For initial data which is expected to disperse, the finite 

dimensional dynamics follow the unstable portion of the soliton curve. 

d 

1. Introduction 

ff^ We consider the initial value problem for the nonlinear Schrodinger equation (NLS) in 
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iut + Au + g{\u\'^)u = 0, 
u{x, 0) = uo{x), 



(1.1) 

where the nonlinearity g{s) is a saturated nonlinearity of the form 



> (1.2) g{s) = s'^ '^'^ 
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C3 where 2+ > p > 2+ ^ > ^ > q > for d > 3 and oo > p > 2 + ^ > ^ > q > hr d < 3. 



For |m| large, (1.1) behaves as though it were subcritical while for |m| small, it behaves 
as though it were supercritical. This guarantees both existence of soliton solutions and 
global well-posedness in H^. 

For our purposes, p must be chosen substantially larger than the critical exponent, 
|, in order to allow sufficient regularity when linearizing the equation. For our numerical 
analysis, we work in one spatial dimension, with the specific nonlinearity 
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The equation (1.1) is globally well-posed in fl with the usual norm 



Illl2 ||1|2 ||||2 

ll'"ll//lnL2(|x|2) = W'^Wm + lFllL2{|a;|2)) 

where is the usual Sobolev space with norm 

II Il2 II ||2 I iirv Il2 

W'^Wm = IfIIl2 + II VM||i2 

and L^(|xp) is the weighted Sobolev space with norm 

II ||2 _ II I I l|2 

ll'"llL2(|a;|2) — \\\X\U\\^2. 

This is commonly referred to as the space with finite variance. The global well-posedness 
initial data in U L^(|a;p) follows from the standard well-posedness theory for semilinear 
Schrodinger equations. Additionally, we assume that uq is spherically symmetric, which 
implies u{x, t) is also spherically symmetric for all t > 0. Proofs can be found in numerous 
references including [15] and [43]. 

A soliton solution of ( |1.1 ) is a function u{t,x) of the form 

(1.4) = e*-V^(x), 

where u > and (pu}{.x) is a positive, spherically symmetric, exponentially decaying solution 
of the equation: 

(1.5) - cu0^ + (?(02 = 0. 

For our particular nonlinearity, for any a; > there is a unique solitary wave solution (j)i^{x) 



to ((L5]), see [8] and [29]. 

For large u the solitons are stable, while for small u they are unstable. A precise stability 
criterion identifying stable and unstable regions is provided in [22] and [36], generalizing 
earlier work on stability in [44], [45]. This amounts to examining the relation uj h- > ||(/)a;||^2, 
defining a soliton curve. Where it is increasing(decreasing) as a function of u, the solitons 
are stable (unstable). Several such curves appear in Figure [l} 

As can be seen numerically in Figure [T| the nonlinearity g spawns a soliton of minimal 
mass. Though certain asymptotic methods can be used to describe the increasing nature 
of the curve as — (multiscale methods) and — >■ oo (variational methods), we forego 
an analytic description of the soliton curve and focus on the minimal mass soliton shown 
to exist in the numerical plot. In [16], Comech and Pelinovsky demonstrated that the 
minimal mass soliton possesses a fundamentally nonlinear instability. They accomplished 
this by finding a small perturbation that forces the solution a fixed distance away from the 
minimal mass soliton in finite time. Their technique reduces to studying an ODE modeling 
the perturbation for short times. For appropriate data, the ODE is unstable. 

We conjecture that though the minimal mass soliton may be unstable on short time 
scales, on a longer time scale the solution will ultimately relax to the stable branch of 
the soliton curve. This conjecture is part of a larger conjecture that solutions which do 
not disperse as t — > cxo must eventually converge towards the stable portion of the soliton 
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Figure 1. Plots of the soliton curves {4>{uj) with respect to uj) for a sub- 
critical nonlinearity, critical nonlinearity, supercritical nonlinearity, and the 
saturated nonlinearity (1.3). The curves for the monomial nonlinearities are 
found analytically, while the curves for the saturated nonlinearities are found 
numerically using the method discussed in Section 4.1 All are in d = 1. Here 
u > .001, as the supercritical and saturated cases diverge as a; — > 0. 



curve. For nonlinearities with a specific two power structure, dynamics of this type were 
observed by Pelinovsky, Afanasjev, and Kivshar, who modelled the behavior of solutions 
near a minimal mass soliton by a second order ODE via adiabatic expansion in lj [31]. By 
contrast, our method uses the full dynamical system of modulation parameters to find a 
4-dimensional system of ODEs which is structured to allow for eventual recoupling to the 
continuous spectrum. This conjecture has also been explored numerically by Buslaev and 
Grikurov for two power nonlinearities in [12], where they found that a solution which is 
initially a perturbation of an unstable soliton tends to approach and then oscillate around 
a stable soliton. 

The purpose of this work is to numerically explore this conjecture. Following [16], we 
break the perturbation into the discrete and continuous parts relative to the linearization of 
the Schrodinger operator around the soliton. The discrete portion yields a four dimensional 
system of nonlinear ODEs. We further simplify the system expanding the equations in 
powers of the dependent variables and dropping cubic and higher terms. 

An obstacle in studying these ODEs is that the signs and magnitudes of the coefficients 
are not self-evident, necessitating numerical methods. We compute these numbers, which 
are intimately related to the minimal mass soliton, using the sine spectral method. The 
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use of the sine function for numerically solving differential equations dates to Stenger [39] . 
It has been successfully used in a wide variety of linear and nonlinear, time dependent 
and independent, differential equations, [4,7,10,14,19,25,32,42]. In this work, we first 



numerically solve (1.5) for the soliton as a nonlinear collocation problem. We then use this 
information to compute the generalized kernel of the operator after linearization about the 
soliton. 

With these coefficients in hand, we numerically integrate the ODE system, plotting the 
results. We find that there are two different types of behavior for the finite dimensional 
system, depending on the initial data. If the initial data represents a solution which our 
nonlinear solver indicates can support a soliton, then we find that the solution is oscillatory. 
It is initially attracted to the stable side of the curve, and, over intermediate time scales, 
proceeds to oscillate around the minimal mass soliton. If we initialize with this type of data 
but with the unstable conditions found in [16], the ODEs initially move in the unstable 
direction but quickly reverse, before commencing oscillation. On the other hand, if we 
begin with initial conditions which are expected to disperse as t oo, our data indicate 
that the finite dimensional dynamics push the solution along the unstable soliton curve 
towards the value uj = rather quickly. This solution matches well to the solution for 



(1.1) with corresponding initial data for as long as the mass conservation of the solution 
allows, after which our model continues to follow the unstable soliton curve but the actual 
solution disperses. In [31], the authors observed similar dynamics, with both oscillatory 
and dispersive regimes. 

These ODEs are an approximation valid on a short time interval. This study is the 
beginning of an analysis to show that perturbations of the minimal mass soliton are at- 
tracted to the stable side of the soliton curve. In a forthcoming work we hope to show 
how the continuous-spectrum part of the perturbation interacts with the discrete-spectrum 
perturbation. Based on the work of Soffer and Weinstein, [38] we expect coupling to the 
continuous spectrum to cause radiation damping, which will ultimately cause the solution 
to have damped oscillations and select a soliton on the stable side of the curve. 

This paper is organized as follows. In section 2, we introduce preliminaries and necessary 
definitions. In section 3, we derive the system of ODEs. In section 4, we explain our 
numerical methods for finding the coefficients of the ODEs. In section 5, we show the 
numerical solutions of the ODEs and explain our results. Finally, in section 6 we present 
our conclusions and plans for future work. An appendix contains details of our numerical 
method for computation of the soliton and related coefficients. 
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2. Definitions and Setup 

For data uq G fl L^(|a;p), there are several conserved quantities. Particularly impor- 
tant invariants are: 

Conservation of Mass (or Charge): 

Q{u) = - j \u\'^dx = - I \uo\'^dx. 
Conservation of Energy: 

E{u)= [ \Vu\^dx- [ G{\u\^)dx= I \Vuo\^dx- I G{\uo\'^)dx, 



where 

G{t) = / g{s)ds. 







Detailed proofs of these conservation laws can be easily arrived at by using energy estimates 
or Noether's Theorem, which relates conservation laws to symmetries of an equation. See 
[43] for details. 

With this type of nonlinearity, it is known that soliton solutions to NLS exist and are 



unique. Existence of solitary waves for nonlinearities of the type (1.2) is proved by in [8 



in M} using ODE techniques and in higher dimensions by minimizing the functional 
with respect to the functional 

V{u) = j [G{\u\^) - '^\u\^]dx. 

Then, using a minimizing sequence and Schwarz symmetrization, one infers the existence 
of the nonnegative, spherically symmetric, decreasing soliton solution. Once we know that 
minimizers are radially symmetric, uniqueness can be established via a shooting method, 
showing that the desired soliton occurs at only one initial value, [29]. 

Of great importance is the fact that := Q{(t>uj) and := Ei^cpS) are differentiable 
with respect to uj. This can be determined from the works of Shatah, namely [34], [35]. 



Differentiating (1.5), Q and E all with respect to cu, we have the relation 

Numerics show that if we plot with respect to u for the saturated nonlinearity, the 
soliton curve goes to cxd as u; goes to or cxd and has a global minimum at some = cj* > 0; 
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see Figure [TJ This will be explored in detail in a subsequent numerical work by Marzuola 
[27]. 

We are interested in the stability of these explicit solutions under perturbations of the 
initial data. 

Definition 2.1. The soliton is said to he orbitally stable if, \/e > 0, 3S > such that, 
for any initial data Uq such that \\uo — (pujW < for any t < 0, there is some G M such 
that \\u{x,t) — e^^(f)^{x)\\ < t. 

Definition 2.2. The soliton is said to he asymptotically stable, if, 35 > such that if 
\Wq — 0a;|| < 5, then for large t, 3lj,9 > such that u{x,t) — e*'^*"'"^0^(x) disperses as a 
solution to the corresponding linear prohlem would. 

Variational techniques developed in [44] , [45] and generalized to an abstract setting in [22] 
and [36] tell us that when 6{u) = + ojQ^j is convex, or 5"{uj) > 0, we are guaranteed 
stability under small perturbations, while for 6"{uj) < we are guaranteed that the soliton 
is unstable under small perturbations. For brief reference on this subject, see Chapter 
4 of [43]. For nonlinearities that are twice different iable at the origin and of monomial 
type at infinity (which would include our saturated nonlinearities), asymptotic stability 
has been studied for a finite collection of strongly orbitally stable solitons by Buslaev and 
Perelman [13], Cuccagna [17], and Rodnianski, Soffer and Schlag [33]. 

At a minimum of Q^^, soliton instability is more subtle, because it is due solely to nonlin- 
ear effects. See [16], where this purely nonlinear instability is proved to occur by reducing 
the behavior of the discrete part of the spectrum to an ODE that is unstable for certain 
initial conditions. 

2.1. Linearization about a Soliton. Throughout this section, we use vector notation, ^, 
to represent complex functions. Any function written without vector notation is assumed 

to be real. For example, the complex valued scalar function u + iv will be written . In 
this notation, multiplication by i is represented by the matrix J = ■ We denote 

by (f^ the complex vector ^*g^^ ' where 0^ is the real profile of the soliton with parameter 

u. For simplicity, we suppress the u subscript, writing in place of 0^. 

For later reference, we now explicitly characterize the linearization of NLS about a soliton 
solution. First consider the linear evolution of the perturbation of a soliton via the ansatz: 

(2.1) u = e''^\$M+pix,t)) 

with p = I J . For the purposes of finding the linearized hamiltonian at (pomega we do 

\Plm J 

not need to allow the parameters 6 and uj to modulate, but when we develop our full system 
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of equations in Section |3j parameter modulation will be taken into account. Inserting (2.1) 
into the equation we know that since is a soliton solution we have 

(2.2) J(p), + A(p) -u;p = -g{<P^)p-2g\c^^)<p' (^^^ +0(|pT). 

(This calculation is explained in more detailed at the start of Section 3.) Here we have 
used the following calculation of the nonlinear terms of the perturbation equation: 

((7(10 + pn(0 + P) - 9{\<P?)4>) = (^?(0' + 2#iJe + pL + pL)(0 + PRe + mrn) " 9{<P^)(p) 

= 9'{<P^){p%e + '^<PpRe + P]m){(t> + PRe + iplm) 

(2.3) + ^/(0')(pL + 20Pi?e + pL)'(0 + PRe + ^PI.^) + h.o.t.s. 

The linear terms will be absorbed into the linearized operator JTiuj-, while the quadratic 
terms are handled explicitly; the 0{p'^) terms in the expansion of the equation around (p^j 
will be denoted by N{uj,p) in the sequel. In this work, after expansion in powers of p, we 
drop all terms of order greater than two. 

We are interested in linearizing this equation: 

JH ( ^"^M + h.o.t.s, 



(2.4) 


gfpRe^ 
\PlmJ 


where 




(2.5) 




(2.6) 


L_ = -A 


(2.7) 


L+ = -A 



Pi. 



m 



L_ 
L+ 



Definition 2.3. A Hamiltonian, Ti is called admissible if the following hold: 

1) There are no embedded eigenvalues in the essential spectrum, 

2) The only real eigenvalue in [—uj,uj] is 0, 

3) The values icj are non-resonant . 



Definition 2.4. Let (NLS) be taken with nonlinearity g. We call g admissible if there 
exists a minimal mass soliton, (pmin, for (NLS) and the Hamiltonian, Ti, resulting from 
linearization about 6 



is admissible in terms of Definition 2J[ 



The spectral properties we need for the linearized Hamiltonian equation in order to prove 

However note that it is sometimes 



stability results are precisely those from Definition 2.3 



possible to numerically solve this sort of problem even if Definition 2.3 does not hold; see, 
for example, [12]. Notationally, we refer to Pd as the projection onto the finite dimensional 
discrete spectral subspace of fl L^(|xp) relative to Ti. Similarly, represents 
projection onto the continuous spectral subspace for 7i. 
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In this work, we must simply assume that g is an admissible nonlinearity. However, this 
assumption is justified by the observed dynamics. Great care must be taken in studying 
the spectral properties of a linearized operator; although admissibility is expected to hold 
generically, certain algebraic conditions on the soliton structure itself must factor into the 
analysis, often requiring careful numerical computations. See [18] as an introduction to 
such methods and the difficulties therein. To this end, in the forthcoming work [28], two 
of the authors will look at analytic and computational methods for verifying these spectral 
conditions. 



2.2. The Discrete Spectral Subspace. We approximate perturbations of the minimal 
mass soliton by projecting onto the discrete spectral subspace of the linearized operator. 
We now describe, in detail, the discrete spectral subspace at the minimal mass. 

Let uj* be the value of the soliton parameter at which the minimal mass soliton occurs. 
It is proved in [16] (Lemma 3.8) that the discrete specral subspace of H at cu* has real 

dimension 4. The functions ei = ^ and 62 = (^q^ generalized kernel of Ti. 

at every uj. Clearly, ci is purely imaginary and 62 is real. In addition to ci and 62, at uj* 
there are two more linearly independent elements of D^j, the purely imaginary 63 and the 
purely real 64. 

Applying [16] (Lemma 3.9), 63 and 64 can be extended as continuous functions of cu in 
such a way that 63(0;) is purely imaginary, 64(0;) is purely real. We write 





a{u!) 



and 



64(61;) 



/9(^; 




with a and /3 real-valued functions. The linearized operator, restricted to this subspace, is 



(2.8) 



/ 


1 





\ 








1 














1 


^0 





a{u!) 


0^ 



where a{uj) is a differentiable function that is equal to at u;*. 

Before proceeding to the derivation of the ODEs, it is helpful to make a minor change 
of basis. Our goal is that, in the new basis, |ei, 62, 63, e4|, ^61,63^ = 0, which will make 
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it easier for us to compute the dual basis. Replace 63 by 



9 



63 



63 - 



(61,63) ^ 



61 



61 



a 



To preserve the relationship 63 = JH^^e4, we need to replace 64 by 



64 



64 - 



(61,63) ^ 



61 



-62 



To preserve the relationship JHi^e^ — 62 + a(a;)e4, we replace 62 by 



62 = 1 + 



(61,63) 
Ik-lIP 



62 



To preserve JH^^e2 = ei, we get that ei must be replaced by 



61 = 



, (61,63) ' 
1 + -n^nrr^ ) 61 



61 



With these substitutions, the JH^^ matrix on remains the same and we obtain the 
relationship ^63,61^ = 0. Prom here on we will assume that we are working with this 

— * 

modified basis and simply take ej :— ej for j — 1,2,3, 4. 

— * — * 

We will define to be the dual basis to the revised el within D^. That is, the are 

— * 

defined by e D^^ and 



id, ej) = S, 
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If we make the change of basis described above, then we can compute the C,j as follows. 
Define D = ||e2|p||e4|p — (62,64)^. Then: 



^3 
^4 



rei, 



D 
1 



-62 



(62, 64) 



D 



-64, 



lesir 

(62, 64) 



^ ^ 62 

62 + —^64- 



As with the e/s. 



for j = 1, 3 and 



for j = 2, 4 to distinguish between vectors and their scalar components. 



3. Derivation of the ODEs 

To derive the ODEs we start with a small spherically symmetric perturbation of the 
minimal mass soliton, then project onto the discrete spectral subspace. Here, we closely 
follow [16]. 

We begin with the following ansatz, which allows theta and u to modulate: 
(3.1) u{t) = e(/o + p(t)). 

Recall, we have assumed u to be spherically symmetric, so no other modulation parameters 
occur. Unlike in [16] we do not assume that the rotation variable 6{t) is identically zero, so 
we need to include 6 modulation in our full ansatz. Note that the derivation which follows 
applies for all nonlinearities in any dimension; specialization is required only to get explicit 
numerical results. This model includes the full dynamical system for spherically symmetric 
data and is designed in such a way that coupling to the continuous spectral subspace could 
easily be reintroduced. The authors plan to analyze the effect of that coupling, which is 
expected to be dissipative, in a future work. 



Differentiating (3.1) with respect to t, we get 
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where we represent differentiation with respect to t by ■ and differentiation with respect to 
the sohton parameter uj by ■'. Plugging the above ansatz into the equation and cancelhng 
the phase term yields 

(3.2) - (cu + e){$+p) + uJ$' + Jp+ A0 + Ap + ^(10 + pn(0 + p) = 0. 
Recall that since is a soliton solution, — u;0 + A0 + 5f(|0p)0 = O, yielding 

(3.3) -d$-{uj + d){p) + CoJ^' + Ap + g{\$ + p\'){$ + p) - (7(02)0 = 0. 

We multiply by J, solve for p, and simplify. At the same time, we collect the Ap and 
—ujp terms with the linear portion of g{\<^ + p|2)(0 + p) — 5'(02)0, which yields JH^j as 
defined in (2.5). The remaining terms of the nonlinearity are at least quadratic in p; recall 



that the quadratic terms are described in (2.3) and denoted N{uj^p). 
Defining Pj{t) as the coefficient of ej{t) in p, we have 



P 



PRe 
Pim 



PiCi + p2e2 + psca + p4e4 + Pc- 



Then, the above calculations give us 



(3.4) 



p = JH^p~ 9 



ejp 



UJ 







+ N{uj,p). 



Taking the inner product of (3.4) with each of the as defined in Section 



applying (2.8) yields the following system: 

9-9{^,,Jp) + {i,,N), 



2.2 



and 





(fl,^ = P2 




(6,^ = P3 


(3.5) 


(fa,^ = P4 




(^4,^ = a{ 



From this point forward in our approximation we drop the pc component as a higher order 
error term. Using the product rule, we solve the left hand side for pi and put the extra 
terms from the derivative of the operator that projects onto the discrete spectral subspace 
onto the right hand side. 

We have, as in [16], that 

Pdp = ^ ejpj + ^ elTijPj - LoPdP'dP, 
where we have implicity defined 

and used that 

PdJ^PcP=-d;PdPdP. 
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This gives 

Pi + ^ = P2 - e{iu Jp) + (6, N) + djii^,, P'^p) - TiiP,), 
P2 + = P3 - m, Jp) + (6, N) + ooiiL PdP) - E r2.'0.), 

(3.6) P3 = P4 - His, Jp) + (is, N) + P',p) - r3,P,), 

Pa = a{u;)ps - e{U, Jp) + (^4, N) + ^{{U P'aP) - Y ^^.^^i)- 

There is also couphng to the continuous spectrum through terms such as (^i , Jp) which we 
omit. This can be included in the error term and is not analyzed in our finite dimensional 
system. 

To make the system well-determined, we must introduce two orthogonality conditions. 
The first is (p, 62) = 0, and the second is (p, Ci) = 0. These represent the choice of uj{t) and 
9{t) respectively that minimize the size of p. These yields P2 = P2 = 0, and pi = pi = 0, 
respectively. 

The reduced system is then: 

e - -^(Ci, Jp) + (6, N) + djii^,, P',p) - J2 TyPi), 
uJ^Ps- 0(^2, Jp) + (6, N) + a;((6, PdP) - Y r^.-'^.-)' 

(3.7) P3 = P4 - ^(fs, Jp) + fe, N) + dji{^„ P',p) - J2 ^SjPj): 

P4 = a(^)p3 - e{U, Jp) + (^4, N) + cl;((e4, PdP) - Yl ^^jPj^ 

In [16], the authors further reduce this system to prove there is an initial nonlinear 
instability. (Note that they have a slightly different system because they have assumed 
that ^ = 0.) We are interested in the dynamics on an intermediate time scale; thus, we 
retain quadratically nonlinear terms in our equations. 

Our notation is as follows. First, we have 

{il,Jp) = (ei,P20' + P4/3) 
^P4{^1,(3), 

since P2 = 0. Denote 

Cl4=(6(^*),/9(^*)), 

which is the highest order term and the only one that will figure into our quadratic expan- 
sion. Notice that this is a real inner product of functions that normally do not appear in 
the same component of the complex vectors, because of the J in the equation. 
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Similarly, we have 

{i2,Jp) = (6,-Pi0-P3a) 

since pi = 0. Denote 

C23 = 

which is again the highest order term. 
Then we have 

&,JP) = (6,P20' + P4/3> 
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since P2 — 0. Denote 
Finally, we have 



since pi = 0. Denote by 

C43 = {^4{uj*),a{uj*)). 
We also write Qij for the term Tij{uj) = {^i, e'-) ai uo = 00* . 

Next, consider the terms (^j, P'^p). These terms are the Cj components of P'^p. We have: 

" 4 4 
PaP'aP = Pa P)e, + p)^'^ 

4 4 

= XI Xl^S ' ^k)pkej + PsPd^s + P^Pd^A 

j=l k=3 
4 4 

= ^ ^{^j: ^k)pkej + P3(ri3ei + Tagea) + P4(r24e2 + ^4464) 

j=l k=3 

= iil: 63)^361 + (S, 64)^462 + (S, 63)^363 + (S, 64)^464 
P3(ri3ei + rggCa) + P4(r24e2 + r44e4) 

= ((el 63) + ri3)p3ei + ((g, 64) + r24)p4e2 

+ ((S> 63) + r33)p3e3 + ((^4, 64) + 144)^464. 



and 
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Therefore, the relevant nonzero terms are 

(ei,i^>") = ((f;,e3) + ri3)p3, 

(e2,P>-) = ((g,e4)+r24)p4. 

We denote 

Pi3 = (f;(^*),e3(^*)), 

Note that some cancellation will occur with the Tij terms that appear separately in the 
system of ODEs, leaving only these Pij terms in the finally system. 

Finally, the terms {^i, N{u!, p)) must be computed. We are only interested in the qua- 
dratic terms, which, according to (2.3) are: 

(3.8) 3J(7'(0^)0pL + 2 J/(02)0VL + J9'i4>')Mm + 2g' {ct?)ct>pRePim. 

Recall that, since pi and p2 are 0, the projection onto the discrete-spectrum of pRe is just 
paCs and the projection onto the discrete-spectrum of pim is just ^464. We now have to 
compute the lowest-order terms of 

(ei,iV(u;,p)). 

The multiplier of P3 in N{lj, p)) is 

r^i33 = (6,(3^7'(0')0 + 2/(02)02)e2). 

Similarly, we define 



ni44 = 


(6, 




<\>(?4l. 


n234 = 


(6, 


-2g\ 


0^)06364), 


^333 = 


(6, 




^2)0 + 2/(0^)0^)6^), 


^344 = 


(6, 






nA34 = 


(^4, 




0^)06364). 



Notice that, as in the computation of the Cjj, these are real inner products between real 
functions that normally appear in different components of the complex vectors. 
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Lastly, we need to estimate a{uj). Recall that a{uj*) = 0, and that a{uj) appears in (3.7) 



multiplied by ps, so we are seeking only the linear term, a{uj) ~ aQ{uj — oo*). We calculate: 

2 

With these assumptions, we conclude the following: 
Proposition 3.1. The quadratic approximation for the evolution of a perturbation of the 



minimal mass soliton, (3.3), ignoring coupling to the continuous spectrum, is 



= P3 + C23^P3 + ^234P3P4 + ^^P2APa, 
(3.9) P3 = P4 - C34^P4 + ?^333P3 + '^344P4 + ^P33P3, 

Pa = ao(^ - ^^*)P3 + C43^P3 + '^434P3P4 + ^P44P4- 

In this system we implicitly assume that ^, {uj — uj*), pi and their time derivatives are 
all of the same order. 

4. Numerical Methods 



From here on, we use numerical techniques to analyze solutions to (3.9). We will work 
in one space dimension, with the specific saturated nonlinearity g{s) = as described 
in the introduction. Though we have a complete description of the generalized kernel of 
JTi, including its size and the relation among the elements, nothing is expressible in terms 
of elementary functions. As this kernel determines the coefficients in our ODE system, we 
numerically compute it, permitting us to subsequently integrate the ODEs numerically. 

The sine function, sin(7rx) / (vrx) was used to compute solitary wave solutions when ana- 
lytical expressions were not readily available in [26] and in the forthcoming [37]. It has also 
been used to study time dependent nonlinear wave equations, [4, 7, 14, 32] , and a variety of 
linear and nonlinear boundary value problems, [9,10,19-21,30]. 

We will use the sine function to estimate the coefficients in three steps: 

• Compute a discrete representation of the minimal mass soliton, (p^^* . 

• Compute discrete representations of the generalized kernel of Ti, i.e. the derivatives 
with respect to u. 

• Compute necessary inner products for the coefficients. 



4.1. Sine Discretization. The problem of finding a soliton solution of (1.5) is a nonlinear 
boundary value problem posed on M. We respect this description in our discretization 
by approximating functions with the sine spectral method. This technique is thoroughly 



explained in [25,40-42] and briefly in Appendix A.l In the sine discretization, the problem 
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remains posed on M and the boundary conditions, that the solution vanish at ±00, are 
naturally incorporated. 

Given a function u{x) : M ^ M, m is approximated using a superposition of shifted and 
scaled sine functions: 

N 



(4.1) 



CM,N{u,h){x) 



N 

Mfcsinc 

k=-M 



X Xk 

h 



UkS{k,h){x), 



-M 



where Xk = kh for k = —M, . . . ,N are the nodes and h > 0. There are three parameters in 
this discretization, h, M, and A^, determining the number of and spacing of lattice points. 
This is common to numerical methods posed on unbounded domains; see [11]. 

A useful and important feature of this spectral method is that, when evaluated at a node, 

(4.2) Cm,n{u, h){xk) = Uk. 

Additionally, the convergence is rapid both in practice and theoretically. See Theorem [T] 
in Appendix A.l for a statement on optimal convergence. 

Since the soliton is an even function, we may take = M. We will thus write 

(4.3) Cm{u, h){x) = Cm,m{u, h){x). 

The symmetry implies U-k = Uk for k = —M, . . . M. We take advantage of this constraint 
in our computations. In addition, we slave h to M in accordance with (A. 8). 

To compute a discrete sine approximation of the ground state, we frame the soliton 
equation as a nonlinear collocation problem. Approximating 0(x) as in (4.1), we seek 
coefficients {Rk} such that 

dlCM{(p,h){xk) - \CMi(p,h){xk) 

+ 9{\CM{(j), h){xk)\')CM{<P, h){xk) = 0, for = -M, . . . , M. 

By satisfying ( 4.4[ ), the discrete approximation solves the soliton equation in the strong 
sense at the nodes, also known as collocation points. This is in contrast to a Galerkin 
formulation, which solves the equation in the weak sense. However, for the one dimensional 
under consideration, sinc-Galerkin and sinc-coUocation lead to the same algebraic system. 

(4.4) yields a system of nonlinear algebraic equations. Let (p be the column vector 



associated with the discrete approximation of 



(4.5) 



Cm(0, h){xk) ^ 



/ 0-A/ \ 



V 



J 



Differentiation of a sine approximated function that is evaluated at the collocation points 
corresponds to matrix multiplication: 

(4.6) dlCM{<P,h){xk)^D^^^$. 
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Explicitly, D^^) jg 



(4.7) 



D 



(2) 



S{3,h){x)\ 



ft2 3 J - 

^=^k - ^ 1 -2(~i)fe-^- . / , 

h? (k-j)'^ J T 1^ 



Using 4.2 



5((|Cm(0, /i)(a;fc)|^)CAf(0, /i)(xfe) = Ca/(5'(|0|^)0, h){xk) = 5'(|0A;|^)0fc- 



Thus 



^(|Cm(0, h)ixk)f)CMi<P, h)ixk) ^ gm')$ 



( 9{\4>-m?)4>-m\ 



\ 9{\<Pm\ )(j)M J 



With these relations, the discrete system is 

(4.8) D'-^^$-uj$+g{\$\'^)$=0. 

It is this equation to which we apply a nonlinear solver, subject to an appropriate guess. 
We discuss an important subtlety in Appendix A. 2 



4.2. Computing the Minimal Mass. Now that we have an algorithm for finding a 
discrete representation of a soliton, we seek to find the value of the soliton parameter for 
the one possessing minimal mass, along with the corresponding discretized soliton. The 
sine discretization has the property that the L^(]R) inner product is well approximated by 



if, 9) 



„ M 

/ fgdx ^ ^ hfkgk = h{f- g). 



Thus, the mass of the soliton can be estimated by 



L2 



h\<p\\ 



Recognizing that = 0(co'), we seek to minimize the functional h\(j){ijj)\^ with respect to 
uo. The argument uj for which the minimum occurs will be uo* . To find the minimal mass, 
we take the derivative, getting a discrete representation of the minimal mass orthogonality 
condition: 



(4.9) 



We solve (4.9) to find u* , computing 



2h(t) ■ 0' = 0. 
h^* in the process. 



The value uo* can be obtained by other algorithms. In the one-dimensional case, the 
soliton equation possesses a first integral, permitting the minimal mass to be computed 
by numerical quadrature and minimization; a comparison of our results and this approach 
appears in Appendix A.4[ Though these approaches are quite accurate for the task of 
computing the minimal mass, they are inadequate at computing the generalized kernel. 
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Thus, we seek to solve the problem consistently by finding the minimal mass for a given 
2M + 1 dimensional approximation of the problem. 

4.3. Discretized Generalized Kernel. Formally, at the minimal mass soliton 0^^*, there 
are four functions associated with the kernel satisfying the second order equations: 

(4.10) = 0, 

(4.11) L+{-<p'^,) = <f)^,, 

(4.12) L_a = -<Pi,, 

(4.13) L+P = a. 



These four functions can also be discretized with sine, as in (4.1). The operators, L±, have 
discrete spectral representations: 

(4.14) ^ L+ = -Z}(2) +COI- diag{(7(<^0}, 

(4.15) L_ ^ L_ = -Z}(2) +u;I- dmg{g{$l) - 2g'{$l)^l}. 

Taking u = (f>, we successively solve for 0', a, and p. A singular value decomposition must 
be used to get a since L_ has a non-trivial kernel. 

Furthermore, we compute discrete approximations of the derivatives of (p, (f)', a, and (3 



taken with respect to cu at to*. The relevant operators are formed analogously to (4.14) and 



(4.15). 



4.4. Convergence. Amongst the many calculations made, the most important is of u*, 
the parameter of the minimal mass soliton. We summarize the results in Table [ij We see 
that robustly converges, achieving twelve digits of precision and u* appears to achieve 
eleven digits of precision. These are consistent with the values in Table |4] from Appendix 



A. 4, where they were computed using a different methods. 

For the purposes of our simulations, we believe we have sufficient precision, approximately 
ten significant digits, for the time integration of our system of ODEs. Some data for the 



convergence of the coefficients appearing in (3.9) is given in Appendix A. 3 



5. Numerical Results 



We explore here the dynamics of the finite dimensional system (3.9) and compare with 
solutions for the full nonlinear PDF (1.1) with corresponding initial data. 

To solve (3.9), we use the stiff solver ode^S from Matlab after properly preparing the 



initial data using the soliton finding codes in Section 4.1 
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Table 1. The convergence of the sine discretization to the minimal mass sohton. 



M 

1 v± 




uo* 


20 


3. 


820771417633398 


0.177000229690401 


40 


3. 


821145471868853 


0.177576993694258 


60 


3. 


821148930202135 


0.1775^ 


^7655985074 


80 


3. 


821149018422933 


0.1775^ 


^8043323139 


100 


3. 


821149022493814 


0.1775^ 


^8063805561 


200 


3. 


821149022780618 


0.1775^ 


^8065432740 


300 


3. 


821149022780439 


0.1775^ 


^8065432795 


400 


3. 


821149022780896 


0.1775^ 


^8065433095 


500 


3. 


821149022780275 


0.1775^ 


^8065432928 



5.1. Nonlinear Solver. In order to determine the accuracy of our results, we also use 
a nonlinear solver to approximate the solutions with a perturbed minimal-mass soliton 
as initial data. For this nonlinear solver, we use a finite element scheme in space and a 
Crank-Nicholson scheme in time. This is similar to the method used in [23]. In brief, we 



discretize our (1.1) by method of lines, using finite elements in space and Crank-Nicholson 
for time-stepping. This method is conservative, though it is not energy conserving. A 
similar scheme was implemented without potential in [3], where the blow-up for NLS in 
several dimensions was investigated. 

We require the spatial grid to be large enough to ensure negligible interaction with 
the boundary. As absorbing boundary conditions for cubic NLS currently require high 
frequency limits to apply successfully, we choose simply to carefully ensure that our grid is 
large enough in order for the interactions to be negligible throughout the experiment. For 
the convergence of such methods without potentials see the references in [3], [1] and [2]. 

We select a symmetric region about the origin, [— -R, -R], upon which we place a mesh of 
elements. The standard hat function basis is used in the Galerkin approximation. We 

allow for a finer grid in a neighbourhood of length 1 centered at the origin to better study 

the effects of the soliton interactions. 



In terms of the hat basis the PDE (1.1) becomes: 



(mj, v) + v^)/2 - i{g{\u\^)u, v) = 0, 
u{0,x)=uo, u{t, x) = Cv{t)v , 



where (■, ■) is the standard inner product, w is a basis function and u, Uq are linear 
combinations of the w's. 
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Since the f 's are hat functions, we have a tri diagonal hnear system. Let /ij > be a 
uniform time step, and let 

Un = ^c^{nht)v 

V 

be the approximate solution at the nth time step. Implementing Crank-Nicholson, the 
system becomes: 

{Un+l - Un, V) + iht (((Wn+l + Mn)/2)^. , V^) 

= iht {g{\{un+i + M„)/2p)(u„+i + Un)/2,v) , Mo = y^^ayV. 

V 

By defining 

Vn = (m„+i +M„)/2, 

we have simplified our system to: 

{yn,v) + i^{iyn)x,V^) = i^{\yn\'^yn,v) + (m„,w). 

An iteration method from [3] is now used to solve this nonlinear system of equations. 
Namely, we set. 

We take = Un and perform three iterations in order to obtain an approximate solution. 



For our problem, we have taken (1.1) with the nonlinearity 





u 


6 


1 + 


U 


4 



u. 



Then, the minimal mass soliton occurs at 

cu* = .177588065433. 

5.2. Results. With the numerical schemes outlined above, we then compare our finite 
dimensional model to the numerically integrated solution with appropriate initial data 
cuq = c<j(0), ao = P3(0), Po = ^4(0) and Oq = 9{0) = for simplicity. In Figures [2| |3| |4| we 
take /3o > and vary ao, ujq. Similarly, in Figures [5| [6| [7], we take Pq < and once again 



vary ao, ujq. Note that we are comparing solutions to the ODEs to solutions of (1.1) with 
the correct initial parameters so that the initial profiles are identical. 

In the situation where /3o > 0, the initial data is expected to allow the admission of 
a solution with a soliton component as t increases. The finite dimensional system shows 
that if we initially perturb the system either towards the stable or the unstable side of 
the curve, the system produces immediate oscillations. Specifically, if the dynamics begin 



to diverge, the higher order nonlinear corrections in (3.9) arrest the solution, resulting in 
fairly uniform oscillations about the minimal mass soliton. See Figures [2| |3} |4j As one 
can see, for initial values ujo ^ 00*, we see a good fit for several oscillations of our finite 
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dimensional approximation to the dynamics of the full solution. As expected, this weakens 
as Uq diverges from u* due to the nature of our approximations in Section [3j We conjecture 
that such oscillations about the minimal mass when coupled to the continuous spectrum 
will lead generically to a damped convergence of the solution towards the minimal mass 
soliton on a long time scale. 

For (3q < 0, and other initial parameters sufficiently small, the initial data is below the 
minimal mass and therefore is expected to disperse as t increases. In this regime, we see 
an interesting phenomenon, which is that the motion of our finite dimensional system is a 
march along the unstable soliton curve towards the value u = 0. Clearly, the conservation 



laws of (1.1) forbid this from happening for the nonlinear solution on long time scales, and 
we see divergence of the nonlinear solution from our solution as t increases. However on 
shorter time scales we see a good fit of the full solution to the predicted finite dimensional 
dynamics; see Figures [5l |6l [7| Note that this occurs regardless of whether we initially 
perturb in the stable or unstable direction. Our findings confirm the regimes predicted by 
the dynamical system model in [31]. 

It remains to briefiy comment on the convergence of our numerical methods. The stiff 
ODE solver, ode45, for the finite dimensional system of ODEs is a standard Runge-Kutta 
method with known stong convergence results documented in a number of introductory 
texts on numerical methods. In addition, the finite element solver for the full nonlinear 
problem has well-established analytic convergence results, see [1]. Hence, the solutions for 
the corresponding systems are known to be accurate representations of the actual contin- 
uous solutions. Though we have not fully justified in this work the spectral decomposition 



used to derive (3.9), the fact that the infinite dimensional dynamics are so well approx- 
imated by the finite dimensional system constructed from these spectral assumptions is 
quite good evidence that this approximation is a valid one. However, as mentioned in 
Section 2.1, investigating the validity of spectral assumptions will be an important topic of 
future research. 



6. Conclusions and Discussion of Future Work 
In this work, we have used a sine discretization method to compute the coefficients of 



the dynamical system (3.9), which is valid near the minimal mass soliton for a saturated 
nonlinear Schrodinger equation. We find that the dynamical system is an accurate approx- 
imation to the full nonlinear solution in a neighborhood of the minimal mass. Moreover, 
we see that there are two distinct regimes of the dynamical system. 

The first regime represents oscillation along the soliton curve. The finite dimensional 
oscillations are valid solutions on long time scales in the conservative PDF, hence we may 
observe long time closeness of our finite dimensional approximation to the full solution of 

(ra. 
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Comparison of finite dimensional approx. to full solution torut = .177588, a - -.005, |i - .001 



Comparison ot finite dimensional approx. to fuli solution foriu - .177588, u - .005, [i = .001 



- finite dimensional approximation 

- infinite dimensionai solution 




- finite dimensionai approximation 

- infinite dimensional solution 




10 20 30 40 50 60 70 



Figure 2. A plot of the solution to the system of ODE's as well as the full 



solution to (1.1) derived for solutions near the minimal soliton for P3(0) > 
and p3(0) < 0, p4(0) > for ujo = .177588, TV = 1000. 



Comparison of finite dimensional approx. to full soiution forot = .17, n = -.005, |i - .001 



Comparison of finite dimensionai approx. to fuii solution fortn - .17, u - .005, |i - .001 



- finite dimensional approximation 

- infinite dimensionai solution 




- finite dimensionai approximation 

- infinite dimensional solution 




Figure 3. A plot of the solution to the system of ODE's as well as the full 
solution to ( |1.1 ) derived for solutions near the minimal soliton for P3(0) > 
and p3(0) < 0, p4(0) > for ujq = .17, N = 500. 



The second regime represents a strong forcing in our finite dimensional system towards 
the point u = in finite time. As ^ 0, the soliton profile becomes small and broad, 
which is essentially indistinguishable from dispersion. Hence, it is fitting that we observe 
these finite dimensional dynamics precisely when the full solution is expected to become 
completely dispersive. The strong forcing regime is only valid on a finite time scale, since 



^2 oo as C(j — > and (1.1) is conservative. Our numerical evidence suggests that 



dispersive dynamics initially move a solution along the unstable portion of the soliton curve, 
until conservation no longer allows such motion. 
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Comparison of finite dimensional approx. to tull solution forio = .1 5, u = -.005, |5 = .001 



Comparison of finite dimensional approx. to full solution ^orw - .15, u - .005, |i - .001 



- finite dimensional approximation 

- infinite dimensional solution 




- finite dimensional approximation 

- infinite dimensional solution 




Figure 4. A plot of the solution to the system of ODE's as well as the 



full solution to (1.1) derived for solutions near the minimal soliton for ao 



P3(0) > and ao = P3(0) < 0, /3o = P4(0) > for ljq = -15, N = 500. 



Comparison ot finite dimensionai approx. to foil soiution foriu^ - .177588. 
1.04 r 




Companson of finite dimensionai approx. to fuii solution form^ - .177588. - .005. - - 
1.05 r 




Figure 5. A plot of the solution to the system of ODE's as well as the 



full solution to (1.1) derived for solutions near the minimal soliton for ao 
P3(0) > and ao = P3(0) <0, po = p^iO) < for uo = .177588, = 1000. 



We cannot numerically verify our conjecture that soliton preserving perturbations of 
unstable solitons dynamically select stable solitons. However, when we begin with pertur- 
bations that are expected to continue to have a soliton component, we see oscillations about 
the minimal mass; this strongly suggests that, through coupling to the continuous spec- 
trum, the oscillations will damp towards a near-minimal-mass stable soliton. This would 
be quite satisfying from a physical perspective as the system would be moving towards the 
configuration of lowest energy in some sense. 
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Comparison of finite dimensional approx. to full solution ^orw - .17, u - -.005, |i = -.001 



Comparison of finite dimensional approx. to full solution form = .17, .005, \'> - -.001 



- finite dimensional approximation 

- infinite dimensional solution 




- finite dimensional approximation 

- infinite dimensional solution 



Figure 6. A plot of the solution to the system of ODE's as well as the 



full solution to (1.1) derived for solutions near the minimal soliton for ao 



P3(0) > and ao = P3(0) < 0, /3o = P4(0) < for lvq = -17, N = 500. 



Comparison of finite dimensional approx. fo full soiutlon fortrj - .15, a - -.005, |i = -.001 



Comparison of finite dimensional approx. to fuii solution foriu = .15, = .005, |^ = -.001 




- finite dimensional approximation 

- infinite dimensional soiutlon 




7 8 9 10 



Figure 7. A plot of the solution to the system of ODE's as well as the 
full solution to ( 1.1[ ) derived for solutions near the minimal soliton for ao = 
P3(0) > and ao = P3(0) <0, po = P^O) < for ujq = -15, N = 500. 



In this result, we felt it worthwhile to first understand the underlying finite dimensional 



dynamics of (3.7), even in an asymptotic setting. However, in the future, we hope to 



give analytic descriptions of the dynamics of small perturbations of unstable solitons on 
global or near global time scales by looking at the finite dimensional dynamics coupled 
to the continuous spectrum dynamics. In the oscillatory regime, this should result in a 
damped decay of the oscillations to a near minimal mass soliton. In the strong forcing 
regime, this should provide a mechanism for mass transfer into the purely dispersive part 
of the spectrum. Likely, using current techniques this analysis can only be truly done 



A SYSTEM OF ODES FOR A PERTURBATION OF A MINIMAL MASS SOLITON 25 

1.8 



1.6 



1.4 

8 
3 

<S 1.2 

TD 
13 

Q- 

E 
< 

1.0 



0.8 



3.9 4.0 4.1 4.2 4.3 4.4 4.5 

Mass, 1101^11^2 

Figure 8. A plot of the maximum amplitude with respect to the norm 
for a saturated nonlinear Schrodinger equation. Computed at M + 1 = 101 
collocation points for u G [0.01, 1.5]. 

in a perturbative setting, though we believe that initial conditions near strongly unstable 
solitons should exhibit similar behavior. Hopefully more powerful techniques will eventually 
be developed for the global study of the stable soliton curve as an attractor of the full 
nonlinear dynamics. 
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Appendix A. Details of Numerical Methods 

A.l. Sine Approximation. Here we briefly review sine and its properties. The texts 
[25,41] and tlie articles [6,40,42] provide an excellent overview. As noted, sine collocation 
and Galerkin schemes have been used to solve a variety of partial differential equations. 

Recall the definition of sine, 

if ^ = 



(A.l) sinc(2;) = 

and for any /c e Z, /i > 0, let 

(A.2) S{k,h){x) = sine 



X — kh 
h 



The sine function can be used to exactly represent functions in the Paley- Wiener class. 
We spectrally represent functions with sine in a weaker function space. First, we define a 
strip in the complex plane, 

(A.3) Dd = {z eC\\lniz\<d} . 
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Then we define the function space: 

Definition A.l. B^{Dd) is the set of analytic functions on Dd satisfying: 

(A.4a) + L^{-d,d) ~ ^(I^D; as t ^ ±00, with a G [0, 1), 

(A.4b) lim \\f{. + iy)\\^^+ hm ||/(- - iy) ||^, < 00. 

Then, we have the following 

Theorem 1. (Theorem 2.16 of [25]) 

Assume f G B^{Dd), p = 1 or 2, and f satisfies the decay estimate 

(A.5) < Ce-"l^l 

// h is selected such that 

(A.6) h = ^7Td/{aM) < min | nd, tt/v 2 1- , 

then 

\\d^f - d^CMif, /i)|lioo < CM("+i)/2g((-V^&JM))_ 

d identifies a strip in the complex plane, of width 2d, about the real axis in which / is 
analytic. This parameter may not be obvious; others have found d = n /2 sufficient. 

For the NLS equation of order 2a + 1, 

IT 



d 



NLS 



^/2uja 

Saturated NLS "interpolates" between second and seventh order NLS. We thus reason it 



is fair to take d = 7r/\/6uj. Though we do not prove that the soliton and the associated 
elements of the kernel lie in these B^[Dd) spaces or satisfy the hypotheses of Theorem [l| 
we use (A.6) to guide our selection of an optimal h. Since the soliton has a = i/tJ, we 
reason that it should be acceptable to take 



(A.7) - ' ^ 



aM V 6ujM 



(A.7) is dependent on both M and u. Were we to use (A.7) as is, it would comphcate 
approximating, amongst other things, the derivative with respect to 00 of the soliton. To 
avoid this, we use a priori estimates on u*, given in Appendix A. 4 Since we know that 
cij* ~ .18 < .25, it is sufficient to take 

d = 71^/2/3. 

Likewise, since lo* > .1, we may take 
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Thus, instead of (A.7), we use 



h 



7ri 



M ' 



We conjecture that this is a vahd grid spacing for all G (.1, .25); our computations are 
consistent with this assumption. 



A. 2. Numerical Continuation. As discussed in Section the discrete system approx- 
imating (1.5) is 

(A.9) = ^0 + ^(0)0 = 0. 

The multiplication in g{(f))(f) is performed elementwise. In order to solve this discrete system, 
we need a good starting point for our nonlinear solver. We produce this guess by numerical 
continuation. 

Define the function 



X 



g[x;T) 



=0 IS 



1 + rx^ 

Note that g{x, 0) is 7th order NLS and g{x, 1) = g{x), saturated NLS. We now solve 
(A.IO) G(0;r) = D(2)0-A0 + ^(0;r)0 = O. 

At r = 0, the analytic NLS soliton serves as the initial guess for computing (f)r=o- 4'-, 
then the initial guess for solving (A.IO) at r = At. We iterate in r until we reach r = 1. 
This is numerical continuation in the artificial parameter r, [5]. This process succeeds with 
relatively few steps of Ar; in fact only 0(10) steps are required. 

A. 3. Convergence Data. Tabl e [2] offers some examples of the robust and rapid conver- 
gence seen in the coefficients of (3.9). These values are all computed at the minimal mass 
soliton. Also see Table [U 

A. 4. Comparisons with Quadrature Methods. The soliton equation may be inte- 
grated once to get 

1 



(A.ll) 



^0^ + ^[0^-log(l + 0^)] =0. 



Equation (A.ll) yields an implicit algebraic expression for the amplitude, 0(0), 

1 



(A. 12) 



4 



[0(0)^ - log (1 + 0(0)^)] =0. 



Using (A.ll) and (A.12), we can express the mass as 

"HO) 



2 
L2 



(x) dx 



(f){x) dx 



P 



A/-^[/-log(l + /)] 



-1/2 



dp. 
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Table 2. The convergence of several coefficients for the ODE system, com- 
puted at uj* . 



M 


ao 


cl4 


pl3 


nl33 


20 


-0.54851448504 


6.04^ 


^29942099 


-6.04^ 


^29942099 


8.79376331231 


40 


-0.553577138662 


6.12521927361 


-6.12521927361 


8.81625889017 


60 


-0.553555163933 


6.12^ 


^11479039 


-6.12^ 


^11479039 


8.81709463059 


80 


-0.553550441653 


6.12^ 


^27391288 


-6.12^ 


^27391288 


8.81713847275 


100 


-0.553549989603 


6.12^ 


^28576495 


-6.12^ 


^28576495 


8.81714173314 


200 


-0.553549934797 


6.12^ 


^28700415 


-6.12^ 


^28700415 


8.81714206225 


300 


-0.553549934793 


6.12^ 


^28700423 


-6.12^ 


^28700423 


8.81714206227 


400 


-0.553549934795 


6.12^ 


^28700421 


-6.12^ 


^28700421 


8.81714206223 


500 


-0.553549934794 


6.12^ 


^28700423 


-6.12^ 


^28700423 


8.81714206227 



Table 3. Value of the coefficients in (3.9) computed with M = 200 



Coefficient 


Value 


933 


-6.61999411752 


9iA 


-12.4582451458 


Ci4 


6.12828700415 


C23 


1.46358108488 


C34 


4.0422919871 


C43 


0.131304385722 


Pis 


-6.12828700415 


P2A 


-17.9305799071 


P33 


6.61999411752 


Pm 


12.4582451458 


ni33 


8.81714206225 


ni44 


1.84068246508 




1.45559877602 


^^333 


-0.792198288158 


^344 


0.013887281387 


^434 


-0.0822482271619 


ao 


-0.553549934797 



Thus, the mass of the soliton with parameter u is 



(A.13) 



110. 



'aJ|lL2 



P 



^P^-^[/-log(l + P^)] 



-1/2 



dp. 



Equations (A. 12) and (A.13) can be used to approximate ui* by numerically minimizing 



( A.13 ). To compute the amplitude of the soliton, we solve (A. 12 ) using Brent's method with 



a tolerance of l.Oe-14. We use the singular integral integrator QAGS from QUADPACK, 
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Table 4. The sohton parameter and mass of the minimal mass soliton com- 
puted by quadrature. Also included is some of the data for the sine method 
appearing in Table [1} 



Algorithm 


LJ* 




fminbound 
Brent 
golden 
sine with M = 100 
sine with M = 200 
sine with M = 400 


0.177588368745261 
0.177587963826864 
0.177587925853761 
0.177588063805561 
0.177588065432740 
0.177588065433095 


3.821149022780204 
3.821149022778472 
3.821149022776717 
3.821149022493814 
3.821149022780618 
3.821149022780896 



which for this problem is, unfortunately, limited to a relative error of 5.0e-12 and an 
absolute error of l.Oe-15. Trying different routines from the optimization module of SciPy, 
[24], we summarize our results in Table |4j There is a spread of 0(le-12) amongst the 
computed minimal masses and a spread of 0(le-7) amongst the u*. These differences are 
consistent with the prescribed relative error of the quadrature, suggesting the precision of 
this approach to computing the minimal mass and associated u is limited by the quadrature 
algorithm. We summarize these computations in Table |4| which also contains data from 
our sine computations. The sine method is consistent with these quadrature methods. 
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